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ABSTRACT 


Analytical equations derived for evaluating linear 
estimators are applied to extended-Kalman filters for 
approximate performance evaluation. Two cases were con- 
sidered, a single known target trajectory and multiple 
target trajectories with given probabilities of occurrence. 
For the multiple-trajectory case, equations are derived for 
the mean and covariance of estimation error in terms of 
the onal. expectations. Two examples are presented 
to compare the use of the analytical equations with Monte- 


Carlo simulation. 
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I. INTRODUCTION 


In state estimation problems one of the tasks is the 
simulation of a filter or estimator after designing it. 
The purpose of the simulation is to determine the perform- 
ance of the filter under actual operating conditions and 
to investigate sensitivity to inaccuracies or approxima- 
tions present in the design assumptions. 

Monte-Carlo simulation is a common simulation algorithm 
which is currently used. An important drawback of the 
Monte-Carlo approach is the large amounts of computer 
time generally required to achieve a reasonable degree of 
accuracy. To reduce the computation time requirements, 
fewer Monte-Carlo runs maybe used with an attendant loss 
of accuracy. Because of these limitations of the Monte- 
Carlo approach, an alternative method, cosisting of a 
set of analytical equations, has been derived and a 
computer algorithm has been established for the evaluation 
of estimation and prediction error statistics for linear 
filters. It is possible to characterize the propagation 
of the means and covariances of estimation error of a 
filter by difference equations. These difference equations 
can easily be solved using relatively small amounts of 
computer time. The only disadvantage of the analytical 
equation approach is that it only applies to linear filters 


with precomputed gain schedules. 
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In practice, one may be confronted by nonlinear 
problems, such as space vehicle re-entry or orbit deter- 
mination problems, or fire-control estimation problems. 
& common approach for nonlinear problems is to use an 
extended Kalman filter and to evaluate its performance 
by Monte-Carlo simulation. 

The main objective of this research is to investigate 
the possibility of approximate evaluation of the perform- 
ance of extended Kalman filters by applying the analytical 
equation approach. In Chapter II of this thesis the ana- 
lytical equations are derived for linear time-invariant 
€stimators. In Chapter III, the application of the ana- 
lytical equations to problems with multiple tracks occurring 
with given probabilities is investigated and an example is 
presented. In Chapter IV, a common filtering approach for 
nonlinear systems is summarized. In Chapter V, the appli- 
cation of the analytical equations to extended Kalman filters 
is discussed for both the single and multiple track cases 
by using a nonlinear example. Chapter VI contains simulation 
results from another example. The simulations were performed 
using both the analytical equations and the Monte-Carlo 
algorithm for comparison. 

The computer programs that were used are given in 


Appendix A. 
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II. ANALYTICAL EQUATIONS 


A. ASSUMPIIONS 
In the development of analytical equations the following 
assumptions have been made [2] : 
meumne true State trajectory X(k) for k=0,1,2,.... 
is known. 
2. The measurement equation is 
2(k)=h(X(k)) + V(k) (25 al) 
where 
2(k) is the q-dimensional measurement vector at 
time t=kT. 
h is a function (which may be nonlinear) of 
X(k). 
V(k) is the measurement noise vector with the 
assumptions 
a: [ vox) | =0 for all k S> 0. 
(The expected value of the measurement 
noise is zero for all time) 
om R(k) for k=j 
(pv). Ef vcxyy(3)? ] = ' for all k,j 20 
for kfj 
(the measurement noise is uncorrelated). 
(c). Since the target trajectory is known 
E[ x(3)W(x)"] =x(3).E[vtk)] =o for all x, j20 
3. The estimator is linear and described by the equation 


Daf )=Rw/e-V+ E(x) [ 206) -CKr/e-1) | (2.2) 


c2 


The prediction equation is 
X(k/k-1)=¢ X(k-1/k-1)+A\ U(k-1) (253) 
where 
X(k/k) is the n-dimensional estimate of X(k) 
. given measurements 

200) ,2(1),...,4(k). 

X(k/k-1) is the predicted value of X(k) given 
Measurements 2(0), 2(1),...,2(k-1). 

g and © are nxn and qxn matrices, respectively, 
and they are assumed to be known. 

U(k-1) is the m-dimensional deterministic forcing 
vector at time t=(k-1)T. 

A is an nxm known matrix. 

G(k) is the nxm gain matrix. 

The gains may be found by any method. For example, 
optimal estimation gains which minimize the sum of the 
variance of estimation error can be obtained by using the 


Kalman equations: 


G(k)=B(k/k-1) CT [ g Blk/e-1) gT+R(x) | 77 (2.4) 
P(k/k)= [1-G(x) ¢ | P(/k-1) (2.5) 
P(k/k-1)= 9 P(k-1/k-1) gp + Q(k -1) (2.6) 


with the initial conditions 

P(0/-1)=Bo =E{ [ x(0)-Xo] +[ X(0)-Xo]” 
K(0/-1)=Xo 

(In this research gains determined by using the 


Kalman equations have been used for simulations.) 
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Where 
Plc /)=ES | ROc/e) - x0) ] + [ KOZ) - x00) ]® ) (2.7) 
which is the covariance matrix of theoretical - 


estimation error, d 


B(K/k-1)=E( [ R(x /e-1)-K(«) ] [KC /e-2) x) ]7) (2.8) 


Q(k) is the covariance matrix of the random forcing 


maput Wik). i.e. 


a {Q(k) for k=j 
B [ w(x) W(5) l= (2.9) 
0 for k/j 
If the estimation equation (2.2) is initialized with 
the value 


X(0/-1) =X. =E X(0) (2.10) 


O 
it can be shown that the optimal estimate X(x/«) is unbiased, 


Dies 
B[ R(k/e) = X(x)] =O (2.41) 
With the assumptions given above one can derive 
difference equations for the mean of estimation error, the 
covariance of estimation error, the mean of N-step prediction 
error and the covariance of N-step prediction error. These 
equations are derived for the time-invariant estimator given 


by equations (2.2) and (2.3) in the following section. 


B. DERIVATION OF ANALYTICAL EQUATIONS 
1. Mean and covariance of estimation error 
The true state of the target at time k is X(k); thus, 
using Equations (2.2) and (2.3) the estimation error X(k) at 


time k is 
14 
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Rk) # X(u/e) - X(c) 
= 2 X(k-1/k-1) +A U(R-1) + CK) E 20) 
- € @ X(k-1/k-1) - CAU(k-1) | (is) |) 2.22) 
Substituting Equation (2.1) into (2.12) gives 
K(x) =[Z - g(x) ¢ ] g Rox-1/e-1) + Ge) V(x) 
+ G(e) h(X(k)) - X(x) 
+ [I - g(x) GJ u(x-2) (2.15) 


By defining the deterministic matrices S(k) and D(k) as 


s(x) =[ - ctx) ¢] g (2.14) 
D(x) =[Z - G(x) G14 U(k-1) + G(x) h (X(e)) 
- X(k) . (2515) 


Equation (2.13) becomes 

K(k) = S(k) X(k-1/k-1) + Dl) + G(k) VU) (2.16) 
The mean of estimation error is defined as 

E{X(x)] = Ef g(x) Roc-a/c-2) + BC) + Gk) V(K)] 
The matrices-S(k) and D(k) are deterministic, hence 

E[ a(x) vix)] = ek) Bf v()] = 0 
and using properties of the expectation operator gives 

E [ X0x) | = g(x) E[ X(x-1/x-1) } + B(x) 
Defining 

flix) = 8 [ F009] 
gives 


Ax) = gx) EL ¥(x-1/e-1) ] + BO) (2.17) 
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Equation (2.16) can also be written as 


K(k) = S(k) X(k-1/k-1) + De) + GC) V(x) 
+ [ S(k) X(k-1) - S() X(x-1)] 


= S(k)[ KCe-1/e-1) - X0c-1)] + BOE) 
+ G(k) V(k) + S(k) X(k-1) 


or 


B(k) = S(k) X(k-1) + D(k) + G(k) V(x) 
+ $(k) X(k-1) (2.18) 


Since X(k-1) is deterministic and 

E [ ¥(x-1)] =A(x-1) 
Equation (2.17) can be written as 

Af (xk) = Sx) A(k-1) + BC) + S(k) X(k-1) (219) 
Equation (2.19) defines the mean of estimation error at time 
k, in terms of the mean of estimation error at time (k-1). 


The covariance of estimation error is defined 


as 
oe H{{ Bah recy. cent 1 
= 8 [ x(x) me) | (220) 
where 
¥(k) = K(x) f(s) (2320) 
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One can obtain a difference equation for Y¥(k) by substi- 
Suis) (2.13), and (2.19). into (2.21), that is, 


XC) = Sc) X(k-1) + oh) + G(k) V(ik)+ S(k)X(k-1) 
mee) L4( 1-1) - 2h) - S(k)X(k-1) 
= 8(k) [ XQc-1) -~A(x-1) ] + oe) vO) 
= Sk) ¥(k-1) + @(k) V(x) (2.22) 
Then 
Bx/n) = E[ ¥(x) ¥(x) | 
= E[ S(k) ¥(k-1) ¥(k-1)" $()"] 
+ E[ S(k) ¥(k-1) V(x)" G(x)" 
+ EB [G(x) Ve) ¥(e-2)7 SC)" ] 
+ E[ G(x) Vk) Wk)" G(x)" 


Using properties of the expectation. operator yields 


B(k/k) = SC) Ef ¥(x-2) ¥Ce-1)"] s(x)? 


+ G(k) R(k) G(x)? (2,23) 
To obtain expressions for E [ ¥(%-1) V(«)* | and E[ V(x) ¥(k-1)"] 
one can start by deriving an expression for Y(0). The 


Refination of. Y(0) is 


H] 


¥(0) = K(0) -££(0) 


=F(0/0)e= X(0) = E[ Xo) | 
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= ¥(0/o) - x(0) - E[ &(0/0)- x(0)] 
Since _ E[X(0)] = X(o), 
¥(0) = K(0/o) - B[ K0/o)] (2. 24) 
But 
X(0/o) = &(0/-1) + g(0)[ 2(0) - ¢ £(0/-1) | 
= [2 - (0) ¢] £(0/-2) + G(0) n(x(0)) 


TPG FORA (6) (2.25) 
Substituting (2.25) into (2.24) gives 
¥(0) = [L- g(o) g] X(0/-1) + g(0) n(Xx(0)) 
+ g(o) v(o) - E¢ [zr -G(o) ¢] X(0/-2) 
+ G0) B(X(0)) + G(o) ¥(0) 
However, 
E[ (0) v(o)] = g(o) E[ v(0)] = 0 
E[G(o) n(x(o)) } = G0) n(x(0)) 
and using properties of the expectation operator gives 


x(o) = [2 - g(o) g] £(0/-2) + go) HCK(0)) 


+ g(o) v(o) - [2 - (0) g] B[ £0/-1)] 


- G(0) 4{X(0)) 


But &(0/-1) is a deterministic, known quantity, so 
E[ £(0/-1)] = &(0/-1) and 


¥(0) = (0) ¥(0) (2.26) 
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Using Equation (2.22) 


eet). CO) +. G(1), va) 


HT 


(1) G(0) ¥(o) + g(2) ¥(1) (2.27) 
and it can be shown that 
fete ke 4-1 
Y¥(k) = a 1T st) GCj) WCj) 4G) VG) (2228) 
j=0O L i=0 
Equation (2.28) shows that Y(k) is a linear combination 


of V(k) with coefficients which are known constant matrices, 


1G 
rt) =), Le (k) V(C) (2229)) 
=0 
where 
k- € -1 
Hy eee Y SGk-i) cif) for 1-0,1,2,...,(k-1) 
i=0 
(2230) 


L.(ie) = G(x) 


From Equation (2.29) it is easily seen that 


Ef x(x-1) vox)? ] = B[ v(x) x(x-1)"] = 0 
Because of the assumption that E [ v(x) me = 0 for kj. 
This result reduces Equation (2.23) to 
P(k/k) = §(k) P(k-1/k-1) Seo 
+ G(k) R() Gli)? (2731) 
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Equation (2.31) is the error covariance propogation 
equation which expresses the covariance of estimation 
error at time k in terms of the covariance of estimation 
error at time (k-1). 
22) Initial conditions 
To use Equations (2.19) and (2.31) one needs to know 
initial condition values. 


The mean of estimation error at time k=0 is 


JZ (0) 


E | X(0/o) - x(o)| 


\} 


Pad 
E [ (0/0) | - x(o) 
Using Equation (2.25) this becomes 


# (0) = EM [Z - G(0) G] X(0/-1) + G(0) n(X(o)) 


+@(0) ¥(0)) - X(0) 
But 
E| %(0/-1) | = X(0/-2) 
E[ G(o) n(x(o))] = (0) n(x(o)) 
and E [ (0) v(o)] = 0 
thus, 


Alo) =[ £-9(0) g] R(0/-1) - X00) 
+ (0) h(x(0)) (2.32) 


Equation (2.32) gives the initial condition for Equation (2.19). 
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The covariance of estimation error at time k=0 is 
Blo/o) = B[ x(0) ¥(0)"] (2525) 
Substituting Equation (2.26) into (2.33) gives 


P(0/0) 


E| G(o) v(o) v(o)? g(0)* ] 


Eroy R(O) Glo)” (2.35) 
Bauation (2.35) gives the initial condition for equation (2.31). 


3. Mean and covariance of prediction error 


The one-step prediction is given by 


Rct1/k) = GX(k/e) +AU(x) (2.36) 
and the N-step prediction, based on the estimate X(k/k) is 
N-1 : 
X (s+) = Pi x(x/e) + Yo PB vleti) (2.37) 
1=0 


Defining the deterministic matrix A(N) as 
N-1 ; 
Pane) | gp A Uae) (2.38) 
i=0 


then 
X(ctN/«) = BY R(/) + ACN) (2.39) 


The N-step prediction error is 


XCstN/ic) X (act) = X(t) (2.40) 
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The mean of the N-step prediction error is defined as 
AA (thie) = E{ (cen /e)] 
= B[ oY Rc/c) + ACN) - XCD] 
=f B| £x/n)] + ACN) - XGeHN) (2.42) 


But 
Acc) = Bf R0s/c) - XC) | 
= E [ X(x/e)] - X(x) 
and 
E[ R(x/)] = Ale) + x(x) (2.42) 


Substituting (2.42) into (2.41) gives 
fi Ocrn/c) = g fic) + gh Xe) + ACN) - XQ) (2.43) 
The covariance of N-step prediction error is 


P(kt+N/k) = 2{ [tenn - Borni/r) | ‘| Ect) 


-Apwuro (2.44) 


Using (2.40) and (2.43) yields 


~ — N & 
XOctN/k) HFA OstN/X) = g mee) C+ ah) - X(}4N) 


- B(x) - pS x(x) - 2s) 
whe 


= BY kfc) - £00) -408)] 
= oN [X(x/e) 0d] (2.45) 
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which, when substituted into (2.44) gives 


| 


Bowtn/x) = B( PN [ RO/e) - x) ] [£Oc/0) 


ll 
us 
2D} 

= 

z 
| 


(2.46) 
Equations (2.43) and (2.46) give the mean and 
covariance of N-step prediction error based on the mean 
end covariance of estimation error at time k. 
In practice, the measurement equations are often 
dinear, i.e. 
Z(k) = € X(k) + V(x) (2.47) 
Then, in the derivation of the analytical equations, h(X(k)) 
can be replaced with € X(k). This will change Equation 
2.15) to 
pik) =[Z- g(x) g|[ QucG-2) - x00] (2.48) 
and Equation (2.32) to 
fio) = [E- 3(0) 2] [ o/-2) - x60) (2.49) 
4, Summary of key equations 
(1). Mean of estimation error 
fEAX) = §(k) M(x-1) D(k)G+*S(k) X(k=1) 9 (2.19) 
(2). Covariance of estimation error 
Blk/k) = SC) Blk-1/e-1) $(xc)* 
+ Gk) RU) Cx)" (220) 
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g (2.1) 
Dik) = [Z - g(x) g JQu(x-1) 
+ G(k)H(X(k)) = XU) (2.05) 
with initial conditions 
A(o)=(z - g(0) ¢] R(0/-1) - x(0) 
+G(0) h(X(0)) (2,82) 
B(o/o) = g(0) R(o) (o)* (2.35) 
(3). Mean and covariance of N-step prediction error 


Aevenfe) = gO) + gM x0c) + BO) 


EEN) (259) 
P(xtn/x) = gY Box/x) [ gN) * (2.46) 
where 
N-1 
a(n) = Yo gi* 2 a y(ichi) (2.38) 
i=0 


If the system is linear and time-varying the g,A 
and € matrices are time dependent and are denoted by g(x), 
A(k), C(k). Since 


X(k/k-1) = O(k-1) X(k-1/k-1) +4 (k-1) U(k-1) (2.50) 


Z(k) = h(X(k)) + V(x) (2.51) 
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it is only necessary to replace the g, Aand Q@ matrices with 
B(k-1), A (k-1) and C(k) in Equations (2.14), (2.43), (2.15), 
(2.46) and (2.38). 

The analytical equations derived in this chapter 
are applicable to linear filters with precomputed gain 
schedules. 

From Equation (2.31) it is seen that the covariance 
of estimation error is independent of the track, however, the 


mean of estimation error is track dependent. 
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III. APPLICATION OF THE ANALYTICAL EQUATIONS 
TO THE CASE OF MORE THAN ONE TRACK 


The equations derived in the previous chapter are based 
on one.known track. Application of the analytical equations 
using one known track has been studied in Reference / 1_/ 
and results have been tabulated. 

Kieit is desired to evaluate the performance of a filter 
for various tracks, a different track can be used for each 
iteration (ensemble member) in a Monte-Carlo simulation. An 
alternative approach based on the previously derived analytical 
equations is: 

G@))= Apply the analytical equations to caleulate the 
mean and covariance of estimation error for each track indi- 
vidually; these are the conditional means and covariances. 

(2). Calculate the overall mean and covariance using 
the results of (1) and relationships involving the means of 
conditional expectation. 


Assume that there are several tracks to be considered; 
Ch) (C2) (a) (n) 
CIS yl) 5 ace 5 KE) Fa a X(k) 


where each track has a given probability of occurrence 


Bye Ppa Pye Pa ’ 1.@> 
(i) 
p[ x00) = x00) ] =v, 
Do aay teerherpE St (3.1) 
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(4) (4) 
The mean {4 (x) and the covariance P(k/k) can be calculated 
(i) ci) 
for the i'th track X(k). Since X (k) has been used for the 
(i) (i) 

calculation of the £¢ (x), d(x) is the conditional expectation, 
LCs 

(i) (i) 

~ ~ 

fA ln) = | Zoe Ace) : x00 | (3.2) 
(1)) 


Which is the mean of estimation error given the track X(k). 
Thus, conditional expectations can be computed for each of 
the tracks. 

To calculate the mean of estimation error, one needs a 


relationship between the conditional means and the mean. 


A. CONDITIONAL EXPECTATION 


The conditional expectation is defined as / 2_/: 
+O eo 


E [ A/Bep, | ef 2 fa f,(a / Beb,) da, .daj.. .da, 


Where 


A is a continuous random vector. 
B is a discrete random vector. 
f(a / B=.) is the probability density function of the 


random vector A given BB; - 
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The expected value of a continuous random vector A is 


Equation (3.4) can also be written in terms of the joint 


probability density function of continuous random vectors, i.e. 


Pe Fo 2) ~ 
aj] -..f ---f a fa plarb) ° da,....da,.db,...db, 
0° tas 06 


(35) 
From probability theory the joint probability density function 
can be expressed in terms of conditional density function. 


Since the random vector B is discrete (point conditioning), 


fy lard) 
f,(a vi B=b.) $ PT BoD, (3.6) 
J 
Ox 
fyp(2sb) = fy(a / Beb,) - P| Bob, | (3.7) 


Substituting Equation (3.7) into (3.5) and replacing the 
integrals with a summation sign for the random vector B 


(since B is a discrete random vector) gives 


P[B=b, | ssi 
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The term inside the brackets in Equation (3.8) isthe 
conditional expectation given by Equation (3.3). Thus, 
Equation (3.8) takes the final form of 
E[ 4 J= » E[a/ Bp, ]- P[3-b,] (3.9) 
J 
Equation (3.9) can be used to define the mean of estimation 


error in terms of individual means calculated for each track, 


Ao (6 
foo = E[X(x)] 
= (i) (i) 
=i | £00 / £00 =X00 | P[ X(i)=X(x) | (3.10) 
a7 
n (i) 
“ a JA) a (32) 
i=1 
where 


Ci) 
P, = P[ X(x)=x(x) | 
From Equation (2.31) it is seen that the covariance of 


estimation error is independent of the track for linear 


Sevetems.) 20 1% Can be calculated once for all tracks. 


B. A LINEAR PROBLEM 


An example of the application of the analytical equations 


to a simplified fire control estimation problem has been 
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presented in f/1_/. In this section, the same problem is 
considered for three tracks with given probabilities of 
occurrence. The simulation has been done using both the 
analytical equations (program 'EVAL') and a Monte-Carlo 
algorithm. 

ine problem is a part of a simplified anti-aircraft fire 
control estimation problem in one dimension. The filter is 
a Kalman filter which is derived by solving Equations (2.2) 
through (2.11) subject to the following assumptions. The 
model of target motion is a point mass moving with constant 


velocity. The state equations are 


X(k+1) = g X(k) + PW(k) (3,202) 
where 
fy Wh 
g = (3223) 
0 1 
me 
[*= : (3.14) 
f T 


T is the time between measurements and equals 1 second. 


The measurement equation is 


Ics 
—~ 
an 
Ww 

Hy 


Cx) + WO) 


N 
-_—_ 
ea 
4 

i 


[2 o| X(k) + V(k) (3.15) 
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where 
X, (xk) is the range of the target. 
X,(k) is the range rate of the tareet. 
V(k) is the measurement noise. 
w(k) is a random forcing input. 
It is assumed that 
E[v(x)] = 0 (3.16) 
E[ v(x)?] = 625 (m)? 
= R (3.17) 
E[w(x)] = o (3.18) 


E { w(x)? ] = 225 (m/sec*) % 
=" 9 (3% 10) 


The filter initial conditions are 


50 x 102 m 
E[ x(0)]} = (3.20) 
-600 m/sec 
And the initial covariance matrix is 
HO7he <0 
P(0/-1) = (3522) 
Ooo? 


Using this value for P(0/-1) the gain and covariance 


equations were solved to determine the gain schedule Q(k), 
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In the simulations three tracks (true trajectories of the 


target) were used with initial conditions 


(1) 60 x 10? m 
X(0) = 
-600 m/sec (3523) 
(2) | 55x10? m 
X{0) = 


-600 m/sec (3.25) 


C3) 70 & 102 m 
(0) 


X = 
-600 m/sec ue) 
The probabilities of occurrence of the tracks are 
(2) 
P, = P(X(k)=x(k)] = 0.3 (3.26) 
(1) 
ne Pf x(x)=x(i)] = 039 (3227) 
(1) 
P, = P[x(k)=K(k)] = 0.4 (3.28) 


Figures 3.3 through 3.6 show a comparison of the results 
obtained by using the analytical equations and Monte-Carlo 
Simulation. Continuous curves represent Monte-Carlo results 
and triangles represent the results obtained using the 
analytical equations. 10,000 Monte-Carlo runs were used 


(3000 runs for the first track, 3000 runs for the second 
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track and 4000 runs for the third track). From Figures 3.3 
through 3.6 it is seen that the two sets of results are very 


close -- especially the covariances of estimation error. 
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IV. ESTIMATORS FOR NONLINEAR MODELS 


The extended-Kalman filter is a commonly used approach 
for nonlinear estimation problems. In this chapter, the 
application of extended-Kalman filters to nonlinear models 
is summarized. 

Given a continuous and nonlinear estimation problem, 
in order to apply linear discrete filtering theory the 


problem is first linearized and discretized. 


A. DISCRETIZATION 
Consider a nonlinear continuous system of state and 
observation equations given by 
met) = t( X(t), Ult), +) (4.1) 
2%) = n(X(+):) + ¥ (+) (4.2) 
The discrete-time representation of these equations has the 
form 
X(k+1) = a(X(k), U(k),k) + W(x) (4.3) 
zie) = @ (XU) +2V (Ck) (4.4) 
Equation (4.3) can be obtained from (4.1) by using the 


relationship 
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which is the Taylor series expension of equation (4.1). 
W(k) is a random input which may also be used to approximate 


the effect of the truncated higher-order terms. 


B. LINEARIZATION 

If there is available a nominal trajectory X(k) and the 
control U(k) is known, one can linearize equations (4.3) and 
(4.4) by expanding them in a taylor series about the nominal 


trajectory X'(k). For the state equations this gives 


= ' o2 X(k) -X'(k) 
ge ya eeuisiae 
+ W(k) + Higher order terms (H.0.T) (4.6) 


Defining the matrix g(k) as 


oa 
B(K) =S5- 1 x" (ac) UC) JK) sos) 


and truncating the higher-order terms, Equation (4.6) 
reduces to 
X(kt1) = glk) XCc) + aX" (ke) WC) yk) 


- g(k) X'(k) + Wk) (4.8) 
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Since X'(k) is a known quantity, the second and third terms 
of the right hand side of equation (4.8) are known and 
deterministic. Defining the vector U'(k) as 
U'(k) = a( X'(k),U(k),k ) - P(k) X'(k) (4.9) 
Equation (4.8) becomes 
X(ke+1) = O(k) X(k) + U'(k) + W(k) (4.10) 
Applying the same procedure to Equation (4.4), the 


linearized form of the observation equation is 


ee) = (Ck) X(k) + C(X"(k) } = H(k) X*(k) + V(k) (4.112) 
where 
Hye 9 8S 
aX | X'(k) (4.12) 


Defining the vector B(k) as 


B(k) 


Cok) y= Bk) 2? (xk) (4.13) 
yields 


Z(k) 


H} 


Mic) Cc) + Blk) + V Or) (4,14) 
Equations (4.10) and (4.14) represent a linear time- 
varying model, U'(k) and B(k) represent bias terms resulting 
from the linearization process. The analytical equations can 

be applied to this model by replacing the U(k-1) term in 
(2.15) with the UXk-1) term in Equation (2.37) and defining 


new measurement equations given by 
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A*(k) = 2(k) = Bk) 


= HC)! Xe), + V Cc) (4.15) 


C. THE EXTENDED-KAIMAN FILTER 
Consider a nonlinear discrete system of state and 


observation equations given by 


X(k+1) = £(X(k),k ) + W(k) (4.16) 
eer = C(°X(K)”) + Vk) (4,17) 


In these equations f and C are nonlinear functions of the 
state variables X(k), W(k) is a random forcing input and V(k) 
is measurement noise with the usual assumptions (an uncor- 


related, zero mean random processes). 
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In order to apply the linear filter equations, Equations 
(4.16) and (4.17) are expanded about the nominal trajectory 
me.) if it is available. In practice, it is possible to 
have an idea about the target trajectory for some problems, 
like satellite orbit determination problems, in which case 
one may predict the target trajectory to be very close to the 
true trajectory and thus, be able to define a nominal 
trajectory for the purpose of linearization. In such cases 


the gain schedule can be computed before estimation, but in 
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other kinds of problems, it may be impossible to obtain a 
nominal trajectory. An alternative is to linearize the 
problem at each time point about the best estimates of the 
states currently available. In this case the gains can only 
be calculated for each sample as the estimates are available. 
This approach is called the extended-Kalman filter. 

Using an estimate to evaluate the linearization Equation 


(4.16) gives 
X(kt1) = B(k) XC) + W(x) (4,18) 
where 


Af 
Me). (4.19) 
g OX 


 (k/x) 


Similarly Equation (4.17) yields 


Z2(k) 


ll 


Hk) X(k) + V(r) (4.20) 


where 


H(k) = 2 _ (4,21) 
OX | X(k/k-1) 
B(x) and H(x) are obtained as indicated in the linearization 
process discussed earlier and evaluated using the estimate as 
the nominal track. The filter estimation equation is 
X(k/k) = X(k/k-1) + G(x) | 2(x) - (X(x/x-1))] (4.22) 
The prediction equation is 


Hc /e-1) = 2¢-X(ic-1/e-1), k-1) (4.23) 
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The gains are obtained by using the relationship 
Ge) = P(w/e-1) YOK)" | HO) BO/e-1) Cx)" 
Soak) line 


where 


the theoretical covariance equation is 
Bw/x) = [ZL - ge) HC) | BOc/e-1) 
The covariance propagation equation is 


Pate eal) = oP (51) Ptk-1/k-1) B(x-1)" cre! Gls 


and the initial conditions are 


R(0/-1) =E [x(0) ] 
- X 
—(0) 
P(0/-1) =P 
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V. APPLICATION OF THE ANALYTICAL EQUATIONS TO 

EXTENDED-KAIMAN FILTERS: A RE-ENTRY PROBLEM 

In this chapter, the approximate evaluation of the 
performance of an extended-Kalman filter using the analyti- 
cal equations is discussed. Since the analytical equations 
are based on linear estimators with pre-computed gain 
schedules, application of these equations to a problem 
with gains which are not pre-computed is an approximation. 

In an extended-Kalman filter the gains are evaluated on- 
line using the estimates as they are produced. This means 
that the gains will vary from run-to-run even if the track 
is the same. This variation is caused by the differences 
in measurement noise which in turn effect the estimates. To 
use the analytical equation approach to evaluate an extended- 
Kalman filter it is assumed that the linearization which 
results in Equations (4.10) and (4.15) is performed using 
the true Aas ge This leads to a pre-computed gain schedule 
that is used to approximate the extended-Kalman filter gain 


schedule. 


A. THE RE-ENTRY PROBLEM 


In order to compare the results of using the analytical 
equations and Monte-Carlo simulation, a particular problem 
was selected which contains significant nonlinearities in 
both the state and observation equations. This problem 


deals with estimation of the aititude, velocity and constant 
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metistic coefficient of a vertically falling body. The 
measurements are taken at discrete time intervals of 

1 second by a radar that measures range in the presence 
of (discrete) white gaussian noise. The geometry of the 
problem is illustrated in Figure 5.1. It is assumed that 


the body is falling vertically. 


FIG. 5.1. Geometry of Re-entry problem. 
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The following definitions are used: 


X,(t) : Altitude 

X,(+) : Velocity (downward) 

m Mass (constant) 

Ch Drag coefficient (constant) 

A Reference area for drag evaluation (constant) 
Pp Mass density of the atmosphere 

mn : : Radar altitude 

M Horizontal distance 

r(t) : True range 


It is assumed that the effect of gravity is negligible. 


The equations of motion are 


: ChA Pp 
X(t) = - each) C5ack) 
2m 


The air density is approximated by the exponential function 


Osea age eel ah (5.2) 
where 

Y= 5x 10° (5.3) 
Defining 

EAC re atl (5.4) 


which is a constant ballistic parameter, the state equations 


become 
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E(t) = - X(t)? x,(t) e7 7% ) (5.5) 


which are of the form 
X(t) = £( X(t) ) >) ae 


The measurement r(t) (range) is given by 
r(t) = ¥Mo+ (x,(+) -H )? (ery) 


and is observed at discrete instants of time so that the 
observed sequence is 
Ps. 2 2 
atk) = YM? + (X,(k) -H)? + v(x) (5.8) 
where v(k) is white gaussian noise with zero mean and 
constant variance 


E [ vox? | =n 


In this case, the output nonlinearity has the form 
q 2 2 
CC) yy w+ (tk) = 8) (5a0) 
It is assumed that 
HO 


M 


H 


100,000 ft 
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and the filter initial conditions are 


K(0/-1) =X, 


3x 10° ft 
Mm 


26° 1:0 ft/sec (5.10) 
3.x lo? ft 


The assumed initial covariance matrix is 


10° 0 0 
P(0/-1) =| 0 4x10® ~—o (5.11) 
0 0 10 


and the variance of measurement noise is 


i=) 16ax tho ‘noni (et) (5.12) 


The true initial conditions of the falling body are 


3x102 ft 
5 ly 
50) = | 2x10 ft/sec 
1x1072_ ett 


Applying Equation (4.5) to Equation (5.5) it can be shown 
that the discretized state equations are 
X(k+1) = a( X(k) ) (5.14) 


- ¥ X, (x) 
(5.15) 


1 2 
X (x) = X5 (x) + = X, (k) Xk) e 


X4 ct1) a5 


Meri = Xk) - K(x) * X3(k) e- oe 4 x50) X4(x) 


ve VX Ck) + (i)? K(k) ° ora (C5, 26) 
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K4(k+1) = X4(X) (5-17) 


Binearizing Equations (5.15) - (5.17) it can be shown that 


the 9(k) matrix is 


Bic) = 92 
ax K(k) 
eee “15 
ai) ee 33 (5.18) 
A A A 


where X'(k) is the nominal trajectory and the eas are 


-_ ' 2 1’ = PAC) 
Be (7K)? x4(n) 0 7 *2 (5.19) 
A = + ye Xe _Yx:! (k) (5.20) 
we 77 p(k) X5(k) e 1 
al 0 
Aj = ee ‘S (x)? dat (k) (5.6211) 
| L 22) tie ayioren! -YX, (k) 
ia a x (x)? x5 (i) VX (xk) x oad X,(k) x (c)e% il 
E aye) K5 (1) ° 9-2) X4(k) (5.22) 
rn ny 2) XU) 
ie 1 2 te (k) XA(k) e 1 
pea) PG cage) e ya 
Pek) & ok)" Be a) (5.23) 
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Eee ERE)? Xi) on 2 YX (x) (5.24) 
1 -= 0 (5.25) 
paz? 9 (5.26) 
A. ee (5827) 


The H(k) matrix is 
ac 


= 3 «| XC) 


Xj (ic) -H 
yee + O00 - 4)? (5,28) 
In the simulation of the problem, the true track was generated 
by using the discretized state equations (5.15) through 
(5.17). The linearization is made about the track in the 
application of the analytical equations and U'(k) is defined 


as 


Pies), = fava (ie) 1) Bley XO) (5.29) 
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If the true track is generated by solving the state 


equations then 


X'(k+1) = a( X'(k) sk) (5.30) 
and 

U'(k) = X(t) - Blk) X" (x) | (5.3) 

etki) = X (ic) - g(k-1) X'(x-1) (5:22) 


If one uses a track other than that generated by solving the 
state equations, then Equation (5.29) must be used for 
U'(k). The matrices (kk), H(k) and the vector U(k-1) can 
be used for application of the analytical equations. In the 
Monte-Carlo simulation of the extended-Kalman filter, one 
must replace x (i) with the estimates R(k/k) in Equations 
mee29) through (5.27) and by X(k/k-1) in Equation (5.28) 
( g(x) must be evaluated about &(k/k) and H(k) must be 
evaluated about X(k/k-1) yx 

Breures 5.5 through 5.10 illustrate the results of the 
analytical equations and the Monte-Carlo simulation of the 
extended-Kalman filter. The continuous curves represent 
the Monte-Carlo results (1000 runs) and the triangles 
represent results from the analytical equations. The 
figures show that the analytical equations have predicted 
better performance than that predicted by the Monte-Carlo 
Simulation. Actually, this observation is not generally 


true. The difference between. results of the two methods 
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FIG. 5.2. Altitude history of re-entry vehicle. 
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FIG. 5.8. Variance of altitude estimation error vs. 
time. 
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FIG. 5.9. Variance of velocity estimation error vs. 
time. 


A-SCALE=5.006+00. UNITS INCH. 
Y-SCALE=1.00€+04 UNITS INCH. 


61 


030 


025 


020 


O15 


QL0 


005 


000 


VARIANCE OF BALLISTIC 
PARAMETER ERROR 


6 (£t)7 


TIME, seconds 


ee A A ae 


FIG. 5.10. Variance of ballistic parameter estimation 
error vs. time. 
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depends on the measurement noise, the nature of the nonlinear- 
ities and the filter initial conditions. For the re-entry 
problem the two results are very close. For example, at 
50,000 ft altitude, the mean altitude error difference 
between the two results is 30 feet. 

In order to investigate the noise dependence of the 
difference between two results, the problem has been solved 
by each method for various covariances of the measurement 
noise and the root mean square of the difference has been 
calculated. The root mean square of the difference between 


two results is defined as 


a A ; 1/2 
a ~~ ri 
ae )escp Aeon Dp de ylaed.-~) B00) | (5.33) 
k=1 
N 172 
: ail 2 
on 660) =t4— y [ VAR. (xk) - VAR (x) | 
k=1 
where 
Bye), VAR, (k) are the mean and variance of estimation 
error calculated by the Monte-Carlo 
simulation, 
Ho(k), VAR,(k) are the mean and variance of estimation 


error calculated by the analytical 
equations, 


N is the number of time points, and 
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Aus (I), AVAR os (I) are the rms values of the 
difference between the two 


results. 


The results obtained have been tabulated in Table 5.1. 
From, Table 5.1 it is seen that the root mean square 
differences between the two results for both altitude and 
velocity estimation error increase as the covariance of 
measurement noise increases. But this is not true for 
X., Which is the ballistic coefficient. Table 5.1 indicates 
that estimation of the ballistic coefficient is relatively 
insensitive to the measurement noise level. 

Figures 5.11 through 5.16 represent the results of 
50 Monte-Carlo runs (continuous curves) and the analytical 
equations. It is seen that there is significant difference 


between 1000 Monte-Carlo run-results and 50 Monte-Carlo 


run-results. (See Figure 5.5-5.10). 
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FIG. 5.11. Time history of the mean altitude estimation 
error for 50 Monte-Carlo runs. 
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Next, the problem was solved for a different target 


track with initial conditions given by 
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iO 25 « 10" £t/sec |= b (5.35) 
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fae results are illustrated in Figures 517 throush 5.22. 

One can see that the analytical equations predict worse 
performance than that predicted by the Monte-Carlo simulation. 
For comparison, the CPU times used for both method are 

given below. 
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B. APPLICATION OF THE ANALYTICAL APPROACH TO THE MULTIPLE 

TRACK NONLINEAR CASE 

The results derived in Chapter III can also be applied 
to nonlinear problems. Equation (3.11) can be used to 
calculate the mean of estimation error in terms of the 
conditional means of estimation error calculated for each 
track. The conditional means are computed approximately by 
using the analytical equations. 

Since the matrices g(k) and H(k) are calculated using 
the true target trajectory in the nonlinear case, one can 
see that the covariance of estimation error will be dependent 
on the tracks (see Equation (2.31). Thus, one needs a 
relationship to compute the covariance of estimation error 
in terms of the conditional covariances of estimation error, 
i.e. 

C1). 2) (a) (n) 
P(k/k) = £( P(k/k), P(k/k),...,P(k/k),.., P(k/k) ) 
’ (5.36) 
(i). ee 
where P(k/k) is the conditional avai of estimation 
a 


error corresponding to the i'th track X(k). The conditional 


covariance of estimation error is defined as 


(i) is (i) 
P(k/k) = Covar | Foo x) =x00 (5.37) 
The relationship of the covariance to the conditional 


covariances is derived in the following paragraphs. 
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C. CONDITIONAL COVARIANCE 


Given a discrete random vector B, the expected value of 
the dependent continuous random vector A can be expressed 


in terms of the conditional expectations as 


e [a] =>) 2[4/%b,| © PL a8) (5.38) 
J 
For simplicity Equation (5.38) can be written as 


E[ a= E 


E[ A/B] n (5.39) 


where the subscript B denotes the expectation with respect 
to the random vector B. 


The conditional covariance is defined as 


acaaan 
WwW 
een 
— 


But 
e-=\fa-s{ay]-[a-eiai]?| (5.141) 


which is the (unconditional) covariance matrix. Replacing 
A with 


[a etal, [a-2{ at] T in Equation (5.39) gives #* 


E }[a-z{at] -[a-sai]? 72 Bebel(s.42) 


**Equation (5.39) can be written in the general case as 


Efe(a)} = E{ Bf e(a)/al j where g(a) may be any function. 
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By adding and subtracting E [are] inside the parantheses, 


Equation (5.42) can be written as 


; ie - xfa/a]} + {2fa/e] - efa a) } 


[2 - =fara]) + [faa -{1]) ral 


=E 


A 


on 2i[ G-2feAp - @ - (@ - e[as])? 
+ (i - sfav]) - Gla] - ofa)? +G[aA] - ea] 
(a - efa/s] )" + Gays] - ea) 
- Glass] - efa])"| y 2) , (5.43) 
Using the properties of the expectation operator in the 
inner expectation, the cross terms can be calculated as 
e {@ -2[a] ) » faa] - efah"f- 2{@ -2[a/e)/ 3} 


‘@[a/a] - ¥[a yn 
(5.44) 


(because ef a ] and e{ a/3 | are not functions of A, they are 


deterministic with respect to A for a given value of B. 
This means Ea] e[a/s | e e{a] and Ea E[ 4/8 | /B \ = E,{ 4/8} ). 


Using the properties of the expectation operator yields 
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B\ (a - =[a/s])/2}= =[a/e] - ef efa] } 
= faye] - efays]) 
=o (5.45) 
Thus, the cross terms drop out in Equation (5.43). Since 
(E [4/2] : E[a}) is not a function of A but it is a function 
ef B, then Equation (5.43) reduces to 
21> r| B{ G@-=[a/s])- @ - efae])7 5} 
+ Glas] - efa]) - G[aye] - ef]? (5.46) 


or 


B 

(5.47) 
The inner expectation in the first term of the right hand 
side defines the conditional covariance (see Equation (5.40), 


hence 
B= B{ Bay \ +B 4 Glas] - afa})-efa/e] - efa))", 
(5.48) 


Equation (5.48) gives the required relationship between the 


covariance and conditional covariances, i.e. 
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P(k/k) = z| Covar { 00/20) =KC%0} | 
X(ic) 


(i) (i) | 
cae | Cu) = Mey) (ps) — px) 7 


= X( 
(4) (i) 
2 zB} Bic/e) | +B] (BG) - M09) 
AO) 
(4) 
(BM (ck) - pt (x _ 
X(k) 
which can also be expressed as 
n (i) nN i) = 
P(k/k) =) P(k/k)*p, + > [ te) - fiCx)| 
I, a 
Cy) = 
Cee me (yy > Sp. (5.49) 


where 


iy} (iz) 
Li (k) and P(k/k) are the conditional mean and 
covariance of estimation error for 
it 
the i°th track X(k), 
U(k) and P(k/k) are the overall mean and covariance 
of estimation error and 
P; is the probabiiaty of occurrence o£ 
a 
tiesi4 the track X(k). 
First, using Equation (3.11) one can compute the mean of 
estimation error P(k). Then using Equation (5.49) the 


covariance can be calculated. 
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The re-entry problem was simulated for two tracks 


where 


(1) 
Py = P {X(k)=X(xk) | 


0.5 (5.50) 


(2) 
Po = P {X(k)=X(k)} 


0.5 (Se Suh) 
with 


3 x.10° Tears 


(1) i 
mGO)= 1 2 x 10 ft/sec (5252) 
1x10 1/t 
[ 3.5 x 10> ft. 
(2) mn . 
x(0)= Zak 0) ft/sec (5.59) 


1 x 107 1/ft 


mreures 5.23 through 5.28 illustrate the results. The 
continuous curves represent the Monte-Carlo simulation and the 
triangles represent the results from the analytical equations. 
From the figures one can see that the 1000-member ensemble 
Monte-Carlo simulation predicts better performance. Actually 
the differences between the two results are small compared 

to the values of the state variables. For example, the 


maximum difference between the two results is 150 feet for 
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the altitude error at 300,000 feet altitude (which is the 
initial difference; at later times the two results are 
closer). The CPU times used for these simulations are 


given below. 


CPU Time 
Monte-Carlo simulation (1000 runs).......3 min 50 sec 
Meamte-Carlo simulation (50 -runs)...++.+s.s.e..0 min 17 sec 
DMM Gical CQUALLONS.2.0.ecsescseeeeeees.0 Min 7 Sec 
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VI. A SUBMARINE TRACKING PROBLEM 


This chapter discusses an example with linear state 


equations but nonlinear measurement equations. It is a 


problem in 


which estimates are made of the states of a 


submerged target: heading, velocity, rest frequency, range 


to the target at the closest point of approach (cpa) and 


distance to cpa. One passive sensor (sono bouy) provides 


Measurements of target heading and frequency LOGT The 


states are 


6. 


P 
fe) 


Range to the target at the closest point 
of approach from te eneoe 

X distance to cpa from the sensor (negative 
before cpa, positive after cpa). 


Target speed. 


Target heading. 


: Target rest frequency. 


It is assumed that the target has constant velocity and 


heading. Figure 6.1 illustrates the geometry of the 
problem. The state equations are / 4 _/ 
= + e ° 5 k 
eee ey Ty. 16. ) 
= 7a + a | eee) 
ee ee) + VT + a Ne é. 


Veet) = Vee) + ag( Ky) 


92 


ofl DMT AOAST Sern anes wl 


—— 


iw olomexe me Beekiee rh tetquie 


ine Poemevuneed tasniinoa fae 


OBNse 
a ie . 
3 P wii A 
eee a a 
——— a 
? 
ez 
) 


So 


sotenivee date | 


f ees 
‘ saa 
ak 4 
i + ~~ lt 
ei 4 
~ 
r Af 
i iad I 


ia LY - ( ot 


@, 
TARGET 


SENSOR 


naG 26. ser eR Ss OK filter geometry. 


93 


BED) (ae Ole) Ree | % ) 
Sipe ol Label Ye!) 


with By through Ec the random forcing terms, hence, 


Wee) = VG» Ys Xp Me) (6.2) 
SS) 


The quantities % ; Vy , >. are random changes in heading, 
Ss S O 
velocity and rest frequency. They are assumed to be independent, 


zero mean, piecewise-constant rates of change. They have 


variances defined by 


xe 2 EL *} 


5 - BLY *| 
Ey - EL Ye “| 


The values for the standard deviations were taken from 
typical maneuvering parameters for the target 


Ee = 100° / 1000 sec = 1.74533 x 10° rad/sec 


0, 


Ly ——1O-kts!/ 1000 see = 5.5555 x 107 yds/sec” (653) 
s 

Ls = On tiz/* 1000 sec "=210)..5 x 107? Hz/sec 
) 


With the expressions for the random forcing terms included, 


the state equation become / 4/7 
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where A is +1 for counter-clockwise rotation about the sensor 


and -1 for clockwise rotation about the sensor. A is needed 


to give the correct sign for a given geometry. 


The angle measurement equation is 
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The frequency observation equation is 
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where (k), V p(k) are measurement noises and V_ is the 
velocity of sound in the medium Oi is assumed to be 1640 yds/ 


sec). The excitation matrix Q(k) can be found as / 4 / 


Q(k) = E{W(k) W(x)" } 


zi) TAEy © ite ahs ke, | 0 Kopi TEg © 0 
eae ea 00 as) aie, ks ee sy 
l qt =e ae oe Seon z rg 0 
ee eo i. eect 5 ee 
Q(k) = SYMMETRIC | Ty" | : | ; 
tes || ce 
oat | 
| Z 2 
Boeos 
(6.7) 
The first en terms involve state related terms. 
The state transition matrix is 
1 0 0 ) 0 
0 saad UY 0 0 
ao oO 1.0 0 (6.8) 
0 0 0) ul 0 
0 0 0 0 i 
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and the linearized measurement matrix is 


af 


oR ox OV 


a Oo (6.9) 


pS Kae) ii ke) 
BC) | R'(k)? 
| 
ay ee 
’ 2 9 ' | ' ’ 
E 00 Ceres gl) Re go _ POR) VE) REC) 
F5(k)-V, Boe |) Foe RG? 
! | 
0. | at 0 
Rye MURR ERD WEN ae Stas ered tie ee 
ee Ce ee | f(x) 
= ° - | Q 
Bae) y R (k) , 0 ) Bee) 
| 
(6.10) 


where X'(k) represents the known track if the analytical 

A 
equations are used and the predicted states (X(k/k-1) ) in 
the actual extended-Kalman filter. f'(k) is the predicted 


frequency measurement in the extended-Kalman filter and the 


true frequency measurement in the analytical equations. 


R' (k) is defined as 


R'(k) “VR, (x) + SO (6a) 


and the covariance of measurement noise matrix is 


2 
g, 0 
R= (6.12) 
oe 
0) Of 
7.61543x 1072 0 
0 TiGe xe 16:2 


The filter was simulated using the analytical equations and 
the Monte-Carlo simulation with the initial conditions given 
below. The time between measurements was assumed to be 


100 seconds. 


3.15 x 10° yds. 
-2.61 x 10? yas. 
K(0/-1) = | 3.38 yds/see (6.13) 
5.30213 Rad. 
500.432 Hz. 
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The P(O/-1) matrix is 


1.6245x10" 0 0 0 0 
0 4.057x10° 0 0 0 
P(0/-1)= 0 0 2.853 0 0 (6.14) 
0 0 0 0.596 oO 
0 0 0 0 0.684 


and the true target initial conditions are 


2992.59 Yds. 
-5295.85 Yds. 
X(0) = 4, 504 Te ee (6.15) 
5.06145 Rad. 
500 la ian 


These initial conditions have been obtained from a method 
eivien-in»/ 4/7.) 

The analytical equations predict the mean of estimation 
error to be oscillatory. Monte-Carlo simulation results have 
indicated that the filter was unstable. The estimation 
error becomes unbounded as the target passes through CPA. 

The reason is that the system has two angle measurement 


equations used depending on the sign of X, A (see Equation 


p 
(6,5).). .s2mce the filter has significant estimation error 
up to CPA, the filter "crosses" the CPA at a different 
measurement time than the track does. Thus, at a certain 


instant the filter uses a different angle prediction equation 
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then the actual measurement equation. This results in 

a large residual in the estimation equation and the filter 
begins to diverge. If the filter has estimates with very 
small error before CPA, the estimator and the target would 
Meross’ the CPA in the same sampling interval. This is a 
very difficult task because the geometry of the system 
indicates that any fluctuations in the target heading largely 


affect the R and X 3 
.  -epa cpa 


In order to force the sign of R opalk/k-2) to change in 
the same sampling interval as the track, frequency measure- 
ments were tested to determine whether they were "up doppler" 
or "down doppler", then the sign of the X gpa(K/k-2) was reversed 


accordingly so that X and ee et) have the same sign. 


pa 
The results of this approach are illustrated in Figures 6.2 
through 6.11 for 1000 Monte-Carlo runs shown as the continuous 
curves. The triangles represent the results obtained from 
the analytical equations. 

Comparing the results of the two methods, one can see 
that they are completely different (especially the covariances). 
There are several possible reasons for this, including: 

(1) The extended-Kalman filter uses the estimates and 
predictions for calculation of the 9(k), H(k) and Q(x) 
matrices, whereas the true target track was used in the 
analytical equations. Thus, it is reasonable that the gain 


schedules will differ largely if the estimates are poor. 
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(2) Since the linearized measurement matrix H(k) has 
been used, the analytical equations do not include the time 
sharing angle measurement equations; as far as the analytical 
equations are concerned there is only one angle measurement 
equation, because the two equations yield the same derivative. 

(3) Finally, the analytical equations do not "see" the 
change that has been made in the filter algorithm to force 
the filter to cross the CPA in the same measurement interval 


with the target in the Monte-Carlo simulation. 
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VII. CONCLUSIONS 


It has been shown that the analytical equations can be 
used to evaluate the performance of a linear filter with a 
pre-computed gain schedule. The results obtained from the 
analytical equations are exact and are only closely approxi- 
mated by Monte-Carlo simulations with large ensemble sizes 
and correspondingly large CPU times. 

For the case of the extended-Kalman filter, it is 
possible to obtain analytical results which approximate 
those obtained by Monte-Carlo simulations. How close the 
results are, depends on the performance of the filter. The 
main difference between the two methods is that one uses the 
true track whereas the other uses the estimates and predictions 
as the nominal trajectory to evaluate the matrices g(k), 

H(k) (and Q(k),R(k) in some cases, too). If the filter 
performs poorly and gives diverging estimates, one can 
expect that the matrices Z(k), H(k), Q(k) and R(k) will be 
completely different in the two methods. This will cause 
different results in the two methods. 

The matrix Q(k) effects only the gain schedule. If it is 
a constant matrix or independent of the states, there will be 
no difference between results obtained by the analytical 
equations and Monte-Carlo simulation. The matrix R(k) has 
two kinds of effects. One is similar to that which Q(k) has, 


the other one is that the R(k) matrix represents the covariance 


JekZ 


eo an, 


of measurement noise, and in Chapter V, it has been shown 
that the "larger" R(k) is the, larger the difference between 
results becomes. Brom the results given in Chapter V, it 
is observed that the initial filter states effect the results 
too. Depending on how close the initial filter states are 
to the true initial states, the time at which the filter 
converges will change. 

In summary, one can use the analytical equations to 
approximately predict the performance of an extended-Kalman 
filter. The accuracy of prediction depends on the nature of 


the problem and the filter. 
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APPENDIX A 
COMPUTER PROGRAMS 


A. PROGRAM "EVAL" 

The Program "EVAL" is a computer algorithm which solves 
the analytical equations derived in Chapter II. The algorithm 
has been written for the special case of the analytical 
equations which is developed for the MK. 86 Fire Control 
System as discussed in Reference /1/ . The estimation 


equation has been assumed as 


Ruf) = Re/e-1) + Se) [ MC) 20) + My (ce) BOc-1) 
+ Mo(k) 2(k-2) - © X(k/e-1)] (4.2) 


and the measurement equation as 


2(k) = H X(k) + V(x) (A.2) 


where H is the true measurement matrix is the measurement 
matrix with synthetic measurements. The matrices M(k), 
M,(k), M(x) are measurement weighting matrices which have 
been used in the MK. 86 Fire Control System. 

In the development of this research these matrices are 


defined as 


Mo («) =] (qxq Identity matrix) (A.3) 

= = 4 
M, (xk) = M(x) = 0 (A .4) 
Qg=H | (A.5) 
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Then Equation (A.1) reduces to the simpler form 


R w/e) = K(x/e-2) + GC) [ 2() - ¢ Kx/e-1)] (4.6) 


With the assumptions given above and starting with the 


m@etinition of the estimation error, i.e. 
X(k) = R(«/k) - X(x) (A.7) 


(where X(k) is the true estimation error vector) 


the following analytical results have been derived in / 1_/: 


i § Covariance of estimation ‘error 


P(k/e) = S4(e) Ble-1/k-1) $5)" + So() ROk) Sox)" 
#$,(k) R(k-1) S)(c)” + So(k) R(k-2) So(k)” 
+A,(k) + A, (k)” + An() + AC)” (A.8) 
2. Mean of estimation error 
HM (x) = $4(«) M (x-1) + Ble) + $4() X(k-1) (4.9) 
3. Initializing values 


B(0/o) = $,(0) R(0) $,(0)" (A.10) 


je (0) = [x - G(0) g | X(0/-1) - [z - G(0) M (0 )H] ¥(0) 
Gia) 


4, Mean of N-step prediction error 


L., (n+k/x) = oN B(x) + PY x(x) - X (NK) (A.12) 
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5. Covariance of N-step prediction error 


ay, a - 

P(nex/k) + gp! Bex/x) [ gl") (A.13) 
where 

$(k) = GC) MOC) 


$,(«) = G(k) M(x) 


So(k) = G(x) M,(x) 


8,00 =[E- oo) g] g 
Ai (ic) = §o(k) §,(k-1) R(k-1) $51)" 
Ad(k) = §4(") [$(k-1) $4(k-2) + $)(-1) ] ROe-2) (4)? 


R(x) =[g,() HB - I] XC) + S)(c) HOX(k-1) + S$o(k) H XC) 


With the assumptions given by Equations (A.3) - (A.5), 
Equations (A.8) through (A.13) reduce to the form of 
equations derived in Chapter II. 


The computer program given in / 1_/ has been used after 
making some small changes and adding three subroutines. The 
program listing is given at the end of this Appendix. In 
order to use the program, one must provide subroutine AK 
which computes the linearized state transition matrix £(k), 
subroutine H(k) which calculates the matrix HKK and 
subroutine UK which calculates the vector We (ea), If the 


matrices Q(k) and R(k) are required to be calculated Online), 
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then subroutines QON and RON must also be provided. It is 
also necessary to give the proper values to certain flags 


for required operation. 


B. MONTE-CARLO SIMULATION 


A listing of the program used for Monte-Carlo simulation 
is also given at the end of this Appendix. The program has 
been set up for use with both linear and extended-Kalman 
filters. One can simulate a linear or extended-Kalman 
filter by selecting the proper flags and supplying the proper 
subroutines. For extended-Kalman filter simulation, one 
must provide subroutine MEAS for simulation of measurements, 
subroutine TRACK for generating the target track (if not 
read in), subroutine HKK for calculation of H(k), subroutine 
AK for calculation of @(k), subroutine CK for calculation 
SiG (X(k/k=1) ), subroutine XPRED for calculation of 
X(k+1/k) = £ ( X(k/k) ), and subroutines QON and RON for 


Q(k) and R(k) if they are not read in. 
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